The short— time self diffusion coefficient of a sphere in a 

suspension of rigid rods 

J. Guzowski/'^'^ B. Cichocki,^ E. Wajnryb,^ and G. C. Abade^ 

^ Max-Planck-Institut fur Metallforschung, 
Heisenbergstr. 3, D-70569 Stuttgart, Germany 
^Institut fur Theoretische und Angewandte Physik, 
Universitdt Stuttgart, Pfaffenwaldring 57, D-70569 Stuttgart, Germany 
^Institute of Theoretical Physics, Warsaw University, 
ul. Hoza 69, 00-681 Warsaw, Poland 
^Institute of Fundamental Technological Research, 
ul. Swi§tokrzyska 21, 00-049 Warsaw, Poland 
(Dated: December 2, 2008) 

Abstract 

The short-time self diffusion coefficient of a sphere in a suspension of rigid rods is calculated in 
first order in the rod volume fraction <^. For low rod concentrations the correction to the Einstein 
diffusion constant of the sphere due to the presence of rods is a linear function of cf) with the slope a 
proportional to the equilibrium averaged mobility diminution trace of the sphere interacting with a 
single freely translating and rotating rod. The two-body hydrodynamic interactions are calculated 
using the so-called bead model in which the rod of aspect ratio p is replaced by a stiff linear chain 
of touching spheres. The interactions between spheres are calculated using the multipole method 
with the accuracy controlled by a multipole truncation order and limited only by the computational 
power. A remarkable accuracy is obtained already for the lowest truncation order, which enables 
calculations for very long rods, up to p = 1000. Additionally, the bead model is checked by filling 
the rod with smaller spheres. This procedure shows that for longer rods the basic model provides 
reasonable results varying less than 5% from the model with filling. An analytical expression for 
a as a function of p is derived in the limit of very long rods. We show that in the first order in 
1/logp the correction to the Einstein diffusion constant does not depend on the size of the tracer 
sphere. The higher order corrections depending on the applied model are computed numerically. 
An approximate expression is provided, valid for a wide range of aspect ratios. 



I. INTRODUCTION. 



Binary colloidal dispersions of hard spheres and rigid rods have been of growing interest 
recentlyJ ^ * ^ * ^ * ^ * ^ ! In the last two years Kang et alW^ reported on a series of experiments on 
self diffusion of macromolecular spheres of different sizes in host fd-virus dispersions. In 
the first worl<P^ the authors studied tracer diffusion in a suspension of freely moving rods. 
They derived a theoretical expression for the long-time self diffusion coefficient, based on the 
Smoluchowski equation for a sphere and a rod, neglecting the hydrodynamic interactions. 
The quantitative agreement with the experiment was obtained only for the large spheres, 
while for the small spheres the theoretical values of the diffusion constant turned out to 
be overestimated. This was believed to be due to the growing relevance of the hydrody- 
namic interactions. In another two papers long time diffusion of small spheres through an 
entangled rod network was investigated.'^ In the theoretical approach the hydrodynamic 
interactions were taken into account with the rod network treated as a screening medium 
and the screening length remaining a free parameter. 

The present work is an attempt to provide a more fundamental description of the hydro- 
dynamic interactions in the studied system. We focus on self diffusion of the sphere for short 
times in dilute rod suspensions. The diffusion coefficient can be expanded in the rod volume 
fraction, with the zeroth order term being the Einstein diffusion constant of an isolated 
sphere and the first order term incorporating hydrodynamic interactions of the sphere with 
a single, freely moving rod. We provide a detailed description of these interactions applying 
the so-called bead model, in which the rod is replaced by a linear chain of touching spheres. 

The simplest bead model, in which the beads are treated like point friction sources, has 
been first developed by Kirkwood and RisemarP^ in 1948 and used in calculations of dif- 
fusion coefficients of wormlike chains. Yamakawa et al.^ took into account the finite size of 
beads and used the modified Oseen tensor, obtained from replacing each bead by a shell 
of smaller subunits. Following this methodology, de la Torre et alW provided a numerical 
scheme for the so-called shell model, in which the beads are sufficiently small to reproduce 
the surface of the studied particle. However, a full hydrodynamic description, apart from 
forces and translational motions of the beads, must take into account higher order multi- 
pole components as well as lubrication effects arising for configurations near contact. As an 
example we mention the works of Durlofsky et a/.,'^ who included a number of low-level mul- 
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tipoles as well as lubrication corrections and Ladd, who performed computer simulations 
incorporating multipoles in Cartesian coordinates. The multipole method used here is based 
on the original solution in spherical coordinates proposed by Schmitz and Felderho f^^ * ^^ * ^'^ * ^^ 
and further developed by Cichocki and co-workers^^ and provides very fast numerical con- 
vergence. It has already found many applications, for example in the problem of mobilities of 
spheres conglomerate^^^*^^^*^ and, more recently, transport coefficients in suspensions.'^^'^SEI] 
In the numerical calculations the multipole truncation order L^ax is introduced and the 
accuracy of the results can be easily controlled by changing Lmax- Also the lubrication 
corrections have been incorporated in the scheme, providing highly accurate results for the 
configurations near contact. Cichocki et al.^^ first noted that only the relative motions shall 
contribute to the correction and proposed a method, which separates out the collective mo- 
tions. In the foregoing work a new version of the numerical code is used, appropriate for 
spheres of different radii and thus enabling better imitation of the particle shapes. We study 
two different kinds of bead models: in the non-filled case the rod is represented by an array 
of identical, touching spheres on a straight line, while in the filled case, additional smaller 
spheres fill the gaps in the chain. 

Apart from the numerical calculations, we provide a detailed theoretical description of 
the sphere-rod hydrodynamic interactions. In the frame of the bead model we analyze the 
two-body mobility matrix and derive an expression for the sphere mobility trace valid for 
large interparticle distances. We provide an analytical result for the diffusion coefficient as 
a function of the rod aspect ratio p in the limit of large p. 

In Section II we start with the description of the system and we derive an expression for 
the diffusion coefficient in first order in the rod volume fraction 0. For low rod concentrations 
the correction to the Einstein diffusion constant of the sphere due to the presence of rods 
is a linear function of with the slope a proportional to the equilibrium averaged mobility 
diminution trace of the sphere interacting with a single freely translating and rotating rod. 
In Section III we show how the coefficient a can be calculated numerically for the bead model 
of the rod using the multipole method. In Section IV we present analytical results, valid for 
long rods and based on the theoretical analysis of the two-particle mobility matrix presented 
in Appendix A. In Appendix B we derive an asymptotic form of the sphere mobility for large 
interparticle distances and long rods, which is used to support the numerical calculations 
presented and discussed in Section V. In Section VI we summarize the results. 
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II. SYSTEM DYNAMICS. 



We consider a system consisting of a hard sphere of radius a and identical rigid rods of 
length L and diameter D performing Brownian motion in a fluid of viscosity 1] and volume 
V. The configuration space X of the system is described by the position of the sphere fq, the 
centers of the rods Ri, . . . , Rat and the orientations of the rods ui, . . . , uat. The Brownian 
motion of the particles can be described by the Smoluchowski equation for the probability 
distribution P{X, t) for the configuration X = (fq, Ri, . . . , Rat, Ui, . . . , u^v) at time t: 

^P{X,t) = VP, (1) 
where T) is the Smoluchowski operator given by 



V 



'X 



ksTfiiX) ■ Vx + f^{X) ■ Vx$J , (2) 

where fc^ is the Boltzmann constant, T is the temperature, $ is the interaction potential 
and a shorthand notation has been used, incorporating summation over the particle indeces. 
The differential operator Vx is defined by 



' X 



(Vro,VR„...,VR^,Vu„...,Vu^), (3) 



where Vro, Vr, denote gradients with respect to the position of the sphere and rod i, respec- 
tively, and means the gradient in the spherical coordinates on the unit sphere referring 
to rod i. For an explicit form of this operator we refer to the work of Jones^ The mobility 
matrix i-i{X) relates the forces JF and torques T exerted by the particles on the fluid to their 
translational and rotational velocities U and i7: 




(4) 



where = (Fq, Fi, . . . , Fat) and U = (Uq, Ui, . . . , Uat) are both 3(A^ + 1) dimensional 
vectors while T = (Ti, . . . ,Tn) and f2 = (fii, . . . ,0,^) are 2N dimensional vectors. The 
part of the mobility matrix referring to the sphere /Xgo is given by the equation 

Uo = ■ Fo, (5) 
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where the dependence on the configuration of the whole system X has been indicated. We 
assume that the system is in equihbrium and concentrate on the Brownian motion of the 
sphere. We define the short-time self diffusion coefficient as the time derivative of the mean 
square displacement, in the limit of short times: 

Ds = ~{[^r,{t)f)\t=o. (6) 

where (...) denotes the equilibrium average. The expression (|6]) can be evaluated using 
the Smoluchowski equation ([T| and some properties of the Smoluchowski operator ([2]). The 
derivation can be found in the work of Pusey.!^ The result is: 

D, = hBT(IvtJi%), (7) 
Next, denoting for simplicity =: /Xqq, we apply the cluster expansion 

N N 

Moo(^) = ^00 (ro) + t^oo (ro, ^i) + ^Y1 ^oo (^o, Xi, Xj) + . . . , (8) 

i=l 1=1 j<i 

where Xi = (Rj, Uj) for i = 1, . . . , and /Xog (ro) is the mobility of a single sphere. Assuming 
translational invariance and isotropy of the fluid, we have /^QQ^(ro) = /UqI, where /xq = 1/Co is 
the mobility of an isolated sphere and Co = Svrr/a is the Stokes friction coefficient. Assuming 
low concentration of the rods and using ([T]) and ([s]) we arrive at the expansion of in the 
rod volume fraction (p: 

Ds = Do{l-a<P + o{<f>)), (9) 

where Do = ksTfio is the Einstein diffusion constant, = vN/V with v being the volume 
of the rod and the coefficient a is proportional to the averaged trace of the two-particle 
mobility fi'^^J, 

a = j dRg{R)TT[fiSliR)], (10) 

where the rod is assumed to be centered at the origin, R is the sphere position and g(R,) is 
the two-particle low concentration distribution function given 
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, 0, sphere and rod overlap , , 

9i^)={ (11) 
1, sphere and rod non overlap 

In the next section we show how the sphere mobility ^qo ^^"^i hence, the coefficient a can 
be calculated numerically using the multipole method. 



III. THE MOBILITY MATRIX OF SPHERE IN PRESENCE OF ROD. 

For the calculation of the hydrodynamic interactions we replace the rod of aspect ratio 
p = L/D hj a stiff chain of p identical touching spheres (FIG. [T| model A). We propose 
also an advanced model (FIG. [l| model B) in which the gaps between the spheres are filled 
with rings of additional smaller spheres. In the following we briefly describe the multipole 
method for the hydrodynamic interaction between spheres. 

Imagine that we put A^* spheres at positions (R^, . . . , R^s), moving with prescribed 
translational velocities = (U^, . . . , U%s) and rotational velocities i?'^ = (f2^, . . . , ^%s), in 
a viscous unbounded fluid characterised by some incident velocity field vo(r). We assume 
that the resulting fluid flow v(r) is governed by the Stokes equationsj^SESl 

- Vp + r^V^v = 0, V ■ V = 0, (12) 

where p = p(r) is the pressure field. The resulting forces JF'^ = (F*, . . . , F^s) and torques 
T'^ = (TJ, . . . ,T%s), exerted by the particles on the fluid, can be related to the particle 
velocities by the resistance matrix!^ ^ 





(13) 



where: 



^0 



(vo(RD,...,vo(R^.)), 



(a;o(RD,...,'^o(R^.)), 



(14) 



and 



(15) 
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model A model B 




FIG. 1: Side view of the two different hydrodynamic models of the rod and the tracer sphere (on 
the left). The small spheres in model B form circular rings (only 2 of 9 spheres in each ring are 
shown) . 

The forces JF* and torques are determined by the boundary conditions on the particle 
surfaces. The problem can be solved by introducing force densitie^^^'^ fj(r) induced on 
surfaces of each particle i = 1, . . . , N'^. For the stick boundary conditions the fluid velocity 
on a surface Si of particle i equals 

w,(r) = U,^ + nt X (r - Rt) (16) 
and we obtain the following equation on the force densities: 



TV 

w,(r) = Vo(r) + c?r'T(r - r') • f,(r'), r G 

7 = 1 



z = l,...,N^ (17) 



where T(r) is the fundamental solution of the Stokes equations ( 12 ) called the Oseen tensoi'^Sl 



and equal to 



T(r) = -^(l + rf), (18) 

with r = |r| and r = r/r. Once the force densities are known, the total forces and torques 
can be calculated: 



rfrf,(r), (19) 
t;r(r-R^) X f,(r), (20) 



where i = 1, . . . , N^. Equation (17) can be solved using the multipole expansion in spherical 
coordinates. By projection onto the complete set of multipole function^ vimaii), with 
/ = 1, 2, . . ., m = — /, cr = 0, 1, 2 and i = 1, . . . , N'^, being the solutions of the Stokes 



equations (12), we obtain an infinite set of linear algebraic equations on the force multipoles 



fimaii)- In matrix representation in the multipole space, Eq. (17) becomes 



w(i) - vo(2) = G(u)f (j) + Zo '(Of (^), z = l,...,N' (21) 

where the contribution from the force density located on the particle i and contributions 
from the other particles have been separated. The multipole matrix elements of the single 
particle friction operator Zo(i) and the propagator G{ij) can be found in Refs. 18 and 20. In 
an abbreviated form, including the summation over the particle indices, the formal solution 
reads 

f =(G + Zo^)-^(w-vo), (22) 

In the last step the Cartesian components of the forces and torques are expressed 
by the force multipoles /imo(*) and fimi{i), respectively, whereas the Cartesian components 
of the relative translational velocities U| — vo(Rj) and the relative rotational velocities 
— ujoilli) are expressed by the velocity multipoles fimo(0 and Vimi{i), respectively. In 



other words, the resistance matrix (13) is obtained by a projection of Eq. (22) onto the 
FT-subspace. 

In the numerical applications the multipole series is truncated and only the multipoles 
with / < Lmax are taken into account. The resulting resistance matrix CL^ax converges 
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quickly with L^ax^ provided there are no relative motions between the spheres. Other- 
wise, for the configurations near contact, the lubrication effects must be taken into account. 
This has been done by Cichocki et. al^ who adapted the lubrication correction to the mul- 
tipole scheme. It has been shown that separating out the collective motions and applying 
the lubrication corrections only to the relative motions leads to faster convergence of the 
multipole series .l^^l 

Having obtained the A^^-particle resistance matrix we can group the spheres into two 
rigid assemblies, one being a single sphere and the second composed of remaining spheres 
and forming a stiff linear chain. The corresponding two-body resistance matrix for sphere 
and rod is obtained by a linear transformation according to the rigid-body constraints 
imposed on the spheres in the rod. Finally, the two-body mobility matrix is calculated by 
an inversion. Taking the translational part referring to the sphere, we obtain the desired 
quantity ^qo- following section we derive an analytical expression for the coefficient a 

in the limit of large aspect ratios p. We show that the leading contribution can be obtained 
from the bead model taking into account only the lowest multipole. 



IV. ANALYTICAL RESULTS FOR LONG RODS. 



A. The Oseen approximation. 

In many cases a realistic calculation of hydrodynamic properties of sphere conglomerates 
requires incorporating a large number of multipoles.'^ However, it is known that in case of 
linear chains, built of p spheres, the leading term in p results from the lowest multipole .'^'^ 
The torque and rotations as well as higher multipoles are neglected, so that the spheres are 
treated as point friction sources. The force acting on a chosen bead can be written as a sum 
of the single particle Stokes drag and the contributions from other beads: 



F| = C[U| - vo(RD -J2^^J ■ F^l ^ = l,...,p (23) 



where Q is the Stokes friction coefficient of sphere i and T^- = T(R| — Rp. Equation (23) is 



known in the literature as the Oseen approximation.'^'^a this model the resistance matrix 
reduces to its translational part denoted by Cij, hJ = 1; • • • 
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F| = EC,-U,^ z = l,...,p (24) 

So far we did not make any assumptions for the configuration of the spheres. In the special 
case of a hnear chain, due to the axial symmetry, assuming equal sizes of spheres so that 
Q = ( for all i, the resistance matrix (^^j can be expressed by two scalar functions (pij and 



C,- = C[0.,dd + V^,,(l-dd)], (25) 
where d is the unit vector parallel to the rod axis. If the beads touch each other, we take 



j)d and the Oseen approximation (23) leads to the following equations for 



the scalar functions: 



k=n 



4 ^ Iz - fcl 

k=n _ 



k 



^kj, 



-i^kj, 



(26) 
(27) 



k=—n,kf^i 

where the beads are now labeled from —n to n and 2r2 + 1 = p. For large p the above discrete 
equations can be replaced by an integral equation of the formP 



S[x-y) = f{x, y) + X I dtK{x, t)f{t, y), 



where K{x,t) = < 



for X — t > 5 

x—t — 

for |a; — t| < 5 
7^ for t — X > 5 

t—x — 



X 



n 



n 



n 



and f{x,y) 



n(f)ij for A = I 
mpij for A = I 



(28) 
(29) 

(30) 
(31) 



The solution of this type of equation was found by Riseman and Kirkwood^ by means 
of the Fourier transform. To a detailed discussion of the method we refer also to the review 
of Zwanzig et al.?^ In the limit p cxd we find 
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fi^, y) = 7TY1 — ^(^ -y) + o((iogp)-2), (32) 

where 5{x — y) is the Dirac delta. We can apply the Oseen approximation in its general 
form ( [23] ) also for a system consisting of an additional sphere outside the chain. This brings 
us back to the rod-sphere problem. Having found the friction functions ipij and ipij and 
using the results from Appendix A, we can calculate the coefficient a, defined in Section II, 
analytically in the limit p 00. 



B. Coefficient a in leading order in p. 

We choose the coordinate system such that the rod axis points in the Z direction and 
the sphere position R is described by the distance i? = |R| and the angle 6 between R and 
the Z-axis (FIG. [2]). The mobility matrix of the sphere can be obtained by summing up 
contributions from all possible scattering sequences between the sphere and the beads in the 
rod, starting and ending on the sphere. In the limit p ^ 00 the leading contribution comes 
from the sequence containing two propagators attached to the sphere (FIG. [2]). Then: 

Mto = ~ ■ Cij ■ Toj, (33) 

where C,^j is the effective resistance matrix for the beads of the force- and torque-free rod 
and 

Toi = T(R-di). (34) 
Due to the axial symmetry C,^^ can be written in the following form: 

4. = C[0^,dd + V^,,(l-dd)], (35) 

where 0jj and ipij are the scalar friction functions modified according to the constraint of 
rigid body motion of the rod (see Appendix A), 



(Pij <Pij ^ , ; 
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FIG. 2: Configuration of tlie sphere interacting witli two arbitrary beads i,j in tlie rod. Toi 
denotes tlie Oseen tensor and Cij is the resistance matrix of a freely moving rod. Bold arrows 



denote consecutive terms in Eq. ( 33 ) . 



where the functions (pij and ipij are defined in Eq. (25) and we have skipped the summation 



boundaries. From Eq. (36) and simple symmetry properties it follows that 



ij 0, 



i j 
i j 



(37) 



We note that the above relations are consistent with the fact that the total force and torque 
on the rod vanish. 



The expression (33) can be used in the calculation of the coefficient a using Eq. (10) 
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The Oseen tensors Toi decay like 0{R~^), so that the integrals diverge, but the diverging 
parts cancel under summation. We can avoid dealing with the infinities, when we replace 
the tensors Toi by To, — T(R), which behave like 0{R~^). This can be done as soon as 



relations (37) hold. 

We take the rod length L as the length unit and introduce dimensionless quantities 



R* := R/L, 
T*(R*) :=KT(R), 

t* ■= ^ 

The volume f of a cap-ended cylinder of aspect ratio p and diameter D equals 



(38) 



ttD' 



P 



(39) 



In the limit of large rod lengths L oo, keeping the rod diameter D and the tracer sphere 
radius a constant, the distribution function g(R,*) equals unity on the whole space and the 



asymptotic form of a, according to Eqs. (10), (35), (38) and (39), reads: 



a 



J rfR*Tr([T*, - T*(R*)] ■ Cl, ■ [T*^- - T*(R*)]), 



Performing the integration leads to 



(40) 



a 



Using Eqs. (31), (32) and (36) and replacing the sums by integrals, we arrive at 



(41) 
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« = TT7r7r^[i + o((iogp)-^)], (42) 

30/ioClogp 

where the higher order corrections 0((logj9)^^) contain also terms 0{p^^), which arise from 
approximating sums by integrals. From Eqs. (42) and ^ it follows, that for very thin rods 
the first order correction to the Einstein diffusion constant, equal —Doa(j), is proportional 
to L^/logp and does not depend on the tracer sphere size, as long as the sphere radius a 
is sufficiently small in comparison with the rod length L. More precisely, this result is valid 
in the limit logL ^ logD with D/a kept constant, which means that we must also have 
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logL logo. However, when this condition is not satisfied, the dependence on a appears 



in a correction to Eq. (42), which then must be taken into account. 



In the next two sections we investigate numerically the case of the tracer sphere equal 



to the beads in the rod, 2a = D. Then /ioC = 1 and the prefactor in Eq. (42) reduces to 
17/30, which we compare with the numerical results for very long rods. 

V. NUMERICAL CALCULATIONS AND DISCUSSION OF RESULTS. 



The integration according to Eq. (10) has been performed using the Gaussian quadra- 
ture method. Due to the rotational and refiectional symmetry of the mobility matrix, the 
integration area could be reduced to the first quarter of the XZ-plane. 

To optimize the calculations, the numerical integrals have been performed separately on 
the three sub-areas Ai, i = 1, 2, 3: 



A, = {{X,Z):1< d{{X, Z),S)< Xi}, 

A2 = {{X, Z):X,< d{{X, Z),S)< X2}, 

A3 = {{X,Z):X2<d{{X,Z),S) and X^ + Z^ < {4pf}, 



where all lengths are normalized to the sphere diameter D and d{{X, Z), S) denotes the 
distance between a point (X, Z) and the rod surface S. To optimize the accuracy of the 
quadratures, Xi have been set such that Xj/Xj+i does not exceed 30. For the distances 
TZ = R/D > Ap, according to the considerations in Appendix B, the integral could be 
performed analytically using the expression: 

Tr[^S(7^, 9)] = -^Ail + Ci cos^ 9 + C2 cos" 9), (43) 

where the coefficients A, Ci, C2 were obtained from a fit to the numerical data at 7^ = 7j9. 
FIG. |3] presents the angle dependence of the two-particle mobility trace normalized to its 
value for 9 = 7r/2 together with the numerical fit, as well as the numerical coefficients as 
functions of p. With growing p the values of Ci and C2 approach -6 and 9, in agreement 



with the theoretical predictions (B19). 
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FIG. 3: a) Angle dependence of the two article mobility trace Tr[/XgQ ] normalized to its value for 
9 = 7r/2. Here p = 100, 7^ = 700 and L max — 1- Circles are numerical data points and the solid 
line is a two-parameter fit 1 + Ci cos^ + C2cos^0, see Eq. (43). b) The plots of the number 
coefficients — Ci and C2. Theoretical values are 6 and 9, respectively. Solid lines are a guide to the 
eye. 
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FIG. 4: The numerical values of the coefficient a as a function of the truncation order Lmax for 
p = 10. 

The accuracy of the results has been found to depend mostly on the truncation order 
Lmax- FIG.|4]shows the results for different L^ax for p = 10. It can be seen that for L^ax = 3 
the error becomes smaller than 0.5% and for L^ax = 1 it is still only about 4%. This fast 
convergence can be attributed to the following features of the system. Firstly, the spheres 
in the rod do not exhibit relative motions and the lubrication corrections are not necessary. 
Secondly, as the spheres are put on a straight line, only the interactions vanishing like 
are long-ranged (in the sense that their one-dimensional integral diverges). Accordingly, 
truncation at L^ax = 1 gives results converging to those at L^ax = 3 with the growing 
rod length. This is in contrast to the case of three-dimensional conglomerates of spheres 
studied in Ref. 17, for which also and terms exhibit the long-range character and 
all multipoles up to Lmax = 3 must be incorporated. 

Due to the limitation of the computer memory, only the values of p up to 10^ have been 
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accessible. Then, as long as logp ^ 6 is not a big number, the asymptotic form inEq. (42) 
should not be expected to be very accurate. Moreover, the typical experimental values of p 
are of the order of 10^ or 10^. Hence, for the quantitative comparison with the numerical 
data we have introduced an approximate formula incorporating the logarithmic corrections 
of higher order: 

a = —, (44) 

30 logp — 7(p) 

where 7(p) is a function depending on the shape details of the rod. From the theoretical 



considerations (42) we can suppose that for large p the function should be approximately 
given by a series in {\ogp)~^. From a fit to the numerical data for 20 < p < 1000 we have 
obtained 

7(p)^2.12-^, (45) 
log p 

which reproduces the numerical values for L^^ax = 3 within 2% accuracy for p > 12. 

The numerical values of a as a function of p are presented in FIG. [5] and TABLE |T} For 
p > 20 the case Lmax = 1 differs by less then 2% from the case Lmax = 3. This means 
that for longer rods it is sufficient to incorporate only multipoles with 1 = 1. However, this 
is still more than in the Oseen approximation, because beside the forces also the higher 
multipoles / = l,o" = 1 and / = l,cr = 2 are incorporated. As shown in FIG. Isla) the two 



parameter fit in Eq. (45), provides an excellent agreement with the numerical data. We 



find a deviatiation of about 20% from the asymptotic form in Eq. (42), represented by the 
dotted curve. As mentioned, this is due to the slow convergence of the logarithmic tails. For 
comparison we present a curve with 7 = 0.9, which is the value of a logarithmic correction 
to the intrinsic viscosity of a suspension of cylinders, obtained numerically by Ortega and 
Garcia de la Torr^. 



In FIG. ^6) we compare the expression (45) with the numerical values of 
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logp-— , (46) 

Again we find a good agreement between the fit and the data. 

FIG. [1] presents two different types of bead models applied in the numerical calculations 
and in FIG. |6]and TABLE |T] we compare the results. The optimal filling is provided, when 
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a) 




1000 




1000 



FIG. 5: a) The numerical values of the coefficient a as a function of the rod aspect ratio p for 
the truncation order Lmax = 1 (circles) together with the expression a = 17p/30(logp — j{p)) for 
j{p) = 2.12 — 4.39/ log p obtained from a fit to the numerical data for 20 < p < 1000 (solid line), 
7 = 0.9 for the problem of intrinsic viscosity taken from Ref. 9 (dashed line), 7 = corresponding 
to the Oseen approximation (dotted line), b) The numerical values of logp — 17p/30a together 
with the numerical fit 7(p). We note a very slow convergence to the limiting value 2.12 (marked 
by the horizontal line). 




FIG. 6: The coefficient a as a function of the rod aspect ratio p for the two models of the rod {A 
and B) and the truncation orders Lmax = 1 and Lmax = 3. In model A the rod is replaced by 
a stiff chain of equal spheres and in model B additional smaller spheres fill the gaps in the chain 
(FIG. [T|. The inset presents the same plots for smaller values of p. The lines are guide to the 
eye. The relative difference between the models diminishes with growing p and is less than 5% for 
p > 12 and L^ax = 3. The values converge with growing Lmax, suggesting that incorporating all 
multipoles would lead to nearly the same results for the two models. 

the rings consist of 9 spheres of diameter D/A. Each of the small spheres touches the two 
neighboring big spheres and the imaginary rod surface. Accordingly, we deal with 10 times 
more spheres than in the non-filled case. Because the number of operations grows like A^^ 
with being the total number of spheres, we were limited in this case to aspect ratios 
not exceeding p = 25. The values of a for the filled case diminish with L^ax, opposite to 
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the non-filled case, so that the difference diminishes with growing truncation order. This 
suggests that our results are very close to the exact values for the cylindrical rod. 



VI. CONCLUSION 

We have calculated the short-time self diffusion coefficient of a hard sphere in a suspension 
of rigid rods by applying the bead model for the rod and using the highly accurate multipole 
method for interactions between spheres. The numerical values converge very quickly with 
the truncation order Lmax, giving practically exact results already for Lmax = 1 for long rods. 
This is on the contrary to the case of the drag coefficients of three-dimensional conglomerates 
of spheres calculated in Ref. 17, where all multipoles up to L^ax = 3 must have been taken 
into account. We have shown that the expression for a as a function of the rod aspect ratio 
p can be derived analytically for large p for the simplest model of spheres treated as point 



friction sources (see Eq. (42)). In the limit of large rods, the correction to the Einstein 
diffusion constant has been shown to be independent of the tracer sphere size. By applying 
two different hydrodynamic models for the rod we have checked that only the higher order 
logarithmic corrections, vanishing like (logp)~^, depend on the rod shape details. For very 



long rods an approximate, semi-numerical expression for this correction (see Eq. 45) has 
been found. The quantitative difference between the models has been found to be less than 
5% for p > 12 with the results converging with growing p. This suggests that our results 
provide bounds for the exact value for the cylindrical rod and that already the basic bead 
model (without filling) provides a reasonable approximation. 

So far, we have obtained results for a sphere of a diameter equal to the rod diameter, but 
the calculations can be performed for tracer spheres of arbitrary sizes. Beside the hard-wall- 
like repulsions, also the other forms of interparticle potentials can be easily incorporated in 
the calculations. 
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APPENDIX A: MOBILITY MATRIX FOR SPHERE AND ROD. 

In this appendix we analyze hydrodynamic interactions between a sphere and a rod. In 
the limit of long rods, using the bead model, we estimate the behavior of the friction tensors 
for large p and obtain the mobility matrix ^^qq of the sphere in leading order in p. Denoting 
the sphere and the rod by indices and 1, respectively, we define the corresponding two-body 
resistance matrix in the Cartesian coordinates!^ 







^Aoo 


Aoi 


Boo 


Bol^ 






Fi 




Aio 


An 


Bio 


Bn 




Ui 


To 




Boo 


Boi 


Coo 


Coi 




rio 


\tJ 




VBio 


Bii 


Cio 


CiJ 




[nj 



(Al) 



where Ajj, Bjj, Bjj, Cij, i,j = 1,2 are the two-particle friction tensors. From the Lorentz 
reciprocal theorem follow the symmetry properties: 



A ■ • — A-* 



(A2) 



In the mobility problem we are interested in an inverse relation. 



/Uo\ /aoo aoi boo boA /Fo\ 
Ui _ aio an bio bn Fi 
^0 boo boi Coo Coi To 

\nj \bio bn Cio Cn/ \Ti/ 



(A3) 



where the mobility tensors SLy, bjj, bj^, Cij, i,j 



1, 2 have symmetry properties analogical 



to those for the friction tensors in Eq. (A2). We note that aoo = A*oo desired mobility 

matrix of the sphere in presence of a rod. First we estimate the elements of the resistance 
matrix using the bead model for large p and keep only the leading terms. Like in Section IV, 
we introduce a dimensionless distance R* = R/L and take the limit p oo with R* = const. 



Each of the tensors in Eq. (Al) refers to a particular physical situation. As an example 



we choose Aoi, which describes the case of a translating rod and a fixed sphere. The fluid 
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velocity field, produced by the rod moving with translational velocity Ui, results in the 
force AoiUi, exerted by the sphere on the fluid. In the bead model this force is given by 
a sum over the beads i = 1, . . . ,p of a product CoToiF|, where Tqi is the Oseen tensor and 
is the force on bead i. The latter is again related to the velocities of the beads by the 
resistance matrix ^^j. Taking all velocities equal to Ui, we finally obtain: 

Aoi = Co$^To.-C,., (A4) 

where further reflections between the tracer sphere and the beads can be neglected when 
p oo. To estimate the asymptotic dependence of Aqi on p, we use the axisymmetric form 



(25) of with the friction functions 0jj and ipij. In the limit p oo, according to Eq. 
(32), the functions differ only by a prefactor, so it is sufficient to estimate the double sum 
4>ij / Rqi^ where the factor l/-Roi comes from the Oseen tensor and R^^ = Roi/L. Using 
the integral representation, we obtain 



J2jr.=^ r^a; I" dyt{x;R*,9)f{x,y), (A5) 
where t{x; R*, 9) = R/Roi is a function that parametrically depends on R* and 6 but not on 



p and the function f{x,y), given by Eq. (32), is 0(1/ log p). Hence, we we can write Eq, 



(|A4j) as 

Aoi = r^Ai + o(l/logp), (A6) 
logp 

where ^oi = -^oilR-*) is a matrix independent of p, which we specify at the end of the 
appendix. Similarly, we obtain 



Aoo 
An 
Boi 
Cn 



Col + 



S.0 



Cplogp 



Aoo + o{l/plogp) 



Cp 

logp 

CoP 
\ogp 

Cp^ 

logp 



•^11 + o{p/ logp) 
B + o{p/ logp) 
C + o(p^/logp) 



(A7) 



Remaining tensors of the resistance matrix in Eq. (Al) are of the following orders: 
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Boo 


= 0(l/p^ logp), 


Boi 


= 0{l/p\ogp), 


Bn 


= 0{p/{\ogpy), 


Coo 


= 0(1), 


Cio 


= 0(1/ log p). 



(A8) 



Now let us go back to the mobility problem from Eq. (A3). We can calculate the mobility 
tensor ago by putting Fi = Ti = Tq = in Eq. (Al) and solving with respect to Fq. We 
obtain a direct relation between Fq and Uq, which we invert and keep the leading terms in 
p. Finally, the two particle mobility of the sphere reads 



f^oo = T^iAoo - Al ■ Al' ■ ^01 - ^ ■ ■ B^] + o{l/p\ogp), (A9) 
Cplogp 



where 



Aoo = lim ^ $^[0.,T*, ■ dd ■ T*^. + ^,,T*, • (1 - dd) • T* 



P-+00 p 



A,, = hm ^ 5^[0.,T*, ■ dd + ^,,TS, ■ (1 - dd)] 
A,, = hm ^ V[0.,dd + ^,,{1 - dd)] 



(AlO) 



C = lim - dd) 

p— >oo p-^ ^ 



where is an antisymmetric matrix defined in Eq. (Bll) and Tq^ is the Oseen tensor 



normalized to (p (see Eq. (38)). 



From Eqs. ^A9v and (AlO), after straightforward manipulations, we obtain Eqs. (33), 



(35) and (36) in the main text. 
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APPENDIX B: LIMIT OF LARGE INTERPARTICLE DISTANCES: ROD IN 
LINEAR FIELD 



In this Appendix we derive a far-field form of the sphere mobihty matrix ^qq. Taking 
the hmit p — oo we calculate the coefficients A,Ci,C2 for the two-particle mobility trace 



in Eq. (43) 



If we assume, that the rod acts on the fluid by the force Fi, the torque Ti and the 
symmetric dipole moment Si, the resulting velocity field vi(R) at a position R relative to 
the rod center, where |R| >> L, can be written in the form of the Taylor serie^ 

vi(r) = Fi ■ T(r) + Ti • [^V x T(r)] + Si : [VT(r)]^ + . . . , (Bl) 
where the superscript S denotes the symmetric traceless part and 



Fi 
Ti 
Si 



dsf(r) 



dv 



dsr X f (r), 
ds [r{{r)f, 



(B2) 



dv 



where f(r) is the surface force density on the rod surface dv. As soon as the rod moves 
freely, the lowest non-vanishing force multipole moment is the symmetric dipole moment 
Si- In the following, we evaluate Si in the limit p oo using the bead model. 

Consider a single rod immersed in a linear field characterised by a constant vector uio, 
vorticity loiq and rate of strain eio- Assume the rod to move freely, such that the total force 
and torque vanish. The symmetric force dipole moment of the rod Si is then given by the 
following set of equations: 






ySiy 



^ /ul-ulo^ 



^1 — <-^10 

-eio 



(B3) 



A 
C H 
H My 

where A, C,H and M are the friction tensors of the rod and the remaining tensors are 
zero due to the axial symmetry. According to the Lorentz reciprocal theorem the following 
relations hold: 
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— Ai3a, Ca/3 — Cpa, Haj3^i — H^aP, Ma/Sfiu — M^^a/S- (B4) 

Due to the rotational invariance around the rod axis the general forms of the friction tensors 
can be written with use of the unit vector d and a few scalar coefficients. Following the 
notation of Kim and KarillaP^, they read 



Caf3 — dadfj + Y*-^ {5a(} — dadfi 



where the rod friction coefficients depend on the shape of the body and the operation 

I 1 {a/3) means taking the traceless part, symmetric in the indeces a and (3. Explicit 

forms of the tensors d^''^ d*^^) and d*^^^ read 



(B5) 



r 



r 



I (a/3) 



, '(a/3) '7~7'('^'^) 



(B6) 



(2) 



(a/3) 



2 5a^dpdy 



(a/3) 



, '{a/3) '7~7'('^^) 

LliQ^Lii 1^ (Jj jj^Ujl/ 



Solving (B3) with respect to Si, we obtain 



Si 



[M - H ■ C"^ ■ H] : eio = -M : e 



10) 



(B7) 



where we have introduced a fourth rank tensor M. In the two-body problem we are inter- 
ested in the velocity field produced by the sphere exerting some prescribed force Fq on the 
fluid. The rate of strain around the rod is then given by eio = [VT(r)|r=R]'^ ■ Fq. Accord- 



ingly, using Eqs. (Bl), (B7) and the conditions Fi = and Ti = 0, the mobility matrix of 
the sphere reads 



Moo = /iol - [VT(r)|r=R]^ : M : [VT(r)|r=R] 



(B8) 



Using the explicit forms in Eqs. (B5) and (B6) we get 
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M = X^d(°) + 



H\2 



2YC 



d(i) + z^'M(2). 



(B9) 



From Eqs. (B9) and (B8) the angle dependence of the two-particle mobility trace Tr[/LtQQ^] 
can be shown to be a linear combination of unity, cos^ 9 and cos^ 9. This is due to the form 
of tensors d*^*\ which consist only 2 or 4 vectors d. The distance dependence is 0{R~^) 



because of the two derivatives of the Oseen tensor in Eq. (B8). From these considerations 



follows the general form in Eq. (43), where the values of the coefficients A,Ci,C2 can be 
found in the limit p oo, once the rod friction coefficients X*^, F*^, y^, y*^ and Z'^ are 
calculated. For this purpose we apply the bead model. The individual forces acting on 
the beads contribute to the total moments and, in the Oseen approximation, we obtain 



[AS 

I 5 



i 

i 

Si = D^^L5-F|, 

i 

where the antisymmetric tensor and the symmetric traceless tensor L5 are given by 



(BIO) 



d„,6n 



(Bll) 



Analogously, we can express the velocity multipoles of the rod in terms of the single particle 



velocity multipoles. Then, comparing with Eq. (B3), we obtain the following approximate 
formulae for the friction tensors: 



ij 



(B12) 
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Inserting Eqs. (Bll) and (25) into Eq. (B12) and comparing with Eq. (B5), we obtain 



ij 

ij 



Z^' = 0. 



(B13) 
(B14) 
(B15) 
(B16) 
(B17) 



The coefficient Z'^ in the Oseen approximation equals zero, because it corresponds to a 



velocity field which vanishes along the rod. Inserting the above expressions into Eq. (B9) 
the tensor M reduces to 



(0) 



(B18) 



We note that the absence of the terms proportional to d*^^) and d*^^^ for very long rods has 
a clear physical interpretation. It can be showrpSl that any linear, symmetric and traceless 
tensor E can be written as a sum X]j=o i 2 ^^'^^ • However, only the rate of strain of the 
form d*^*^) : E stretches the rod along its axis not producing a vorticity. Hence, the absence 
of terms proportional to d*^^^ and d*^^) means that effectively, a very long, freely moving rod 
immersed in an arbitrary linear field disturbs the fiuid fiow such, as if it was immersed just 
in a field of a rate of strain d^") : E. 

Here, we indicate a correspondence between the above expression and the general formula 



(33). Namely, the Oseen tensors Toi in Eq. (33), from the summation rules (37), can be 



replaced by Toj — T(R), which for large distances becomes iDd ■ VT. Then, again according 
to Eqs. (37), we are left with /Xqq proportional to in agreement with Eq. (B18). 



(2) 

The two-particle mobility trace Tr[^QQ ] for large distances R can be obtained from Eq. 



(B8). After performing the trace and calculating the double sums in Eq. (BIS) using Eqs. 



(31) and (32), we arrive at 



(2)i 



1287^4 logp' 



6cos^^ + 9cos^0] • [1 + 0((logj9)"^)]. 



(B19) 
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where TZ = R/D. Comparing with (43) we find 



A 



— > 

p^oo 1281ogp 



-6 



P—+00 



Co .9. 



(B20) 
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TABLE I: The first virial coefficient a as a function of rod aspect ratio p for different hydrodynamic 
rod models {A,B) and truncation orders L^ax- 



P a 

A A B B 





■^max — 1 


^max — 3 


L'max — 1 




1 


1.07 


1.83 


1.07 


1.83 


2 


1.33 


1.60 


1.64 


1.68 


3 


1.46 


1.65 






4 


1.63 


1.79 


2.01 


1.88 


6 


2.03 


2.14 


2.44 


2.26 


8 


2.44 


2.54 


2.87 


2.68 


10 


2.84 


2.93 


3.30 


3.09 


12 


3.24 


3.43 


3.71 


3.50 


14 


3.63 


3.72 


4.11 


3.90 


16 


4.01 


4.11 


4.50 


4.29 


20 


4.75 


4.85 


5.26 


5.06 


25 


5.64 


5.75 


6.18 


5.98 


30 


6.48 


6.61 






40 


8.10 


8.25 


8.70 




60 


11.1 


11.3 


11.8 




80 


13.9 


14.1 


14.6 




100 


16.5 


16.8 






120 


19.1 


19.4 






150 


22.7 


23.1 






200 


28.5 








300 


39.3 








500 


59.3 








1000 


104.3 
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